function r = solvemodelnew(sbar,p)
% this function solves for the skill threshold stilde, automation
% thresholds xlow and xup, and labor assignment, wages, output and capital
% share
% for low-skill, high-skill and no automation, the model is solved by
% applying root finding algorithms to equilibrium conditions; for interior
% automation, the model is solved by maximizing net output (we experimented
% with both approaches in each case and found these to be the fastest in
% the respective scenarios)

%% set optimization options for net output maximization
options = optimoptions('fmincon','Display','off','StepTolerance',eps,...
    'Algorithm','interior-point','ConstraintTolerance',p.contol);

%% check for no automation

% solve model without capital
rno = laborassignment(sbar,sbar,0,0,p);

% compute effective cost of labor relative to capital
logomega = log(rno.w) + logcapprod(rno.x,p) - p.q - loglaborprod(rno.sgrid,rno.x,p);

% compute q_0 (where automation starts)
[qzero,indmin] = min(-logomega);

% check whether no automation condition is satisfied (q>q_0) and collect
% output if yes
if p.q < qzero
    
    r.sgrid = rno.sgrid;
    r.sbottom = [];
    r.stop = [];
    r.stilde = 0;
    r.x = rno.x;
    r.xbottom = [];
    r.xtop = [];
    r.xlow = 0;
    r.xup = 0;
    r.w = rno.w;
    r.wbottom = [];
    r.wtop = [];
    r.alphak = rno.alphak;
    r.y = rno.y;
    r.c = rno.c;
    r.qzero = qzero;
    r.exitflag = 1;
    r.exitflag1 = "no automation";
    r.exitflag2 = [];
    
    return

end

%% check for type of automation and solve model with automation

% check for global comparative advantage of least and most skilled workers in most complex
% tasks relative to capital (condition 2)
compadleastskilled = (loglaborprod(0,1,p)-logcapprod(1,p))>(loglaborprod(0,0,p)-logcapprod(0,p));
compadsbar = (loglaborprod(sbar,1,p)-logcapprod(1,p))>(loglaborprod(sbar,0,p)-logcapprod(0,p));
compadmostskilled = (loglaborprod(1,1,p)-logcapprod(1,p))>(loglaborprod(1,0,p)-logcapprod(0,p));

% distinguish cases according to Proposition 3

if (compadsbar == 0) && (compadmostskilled == 1) %interior automation
    
    % compute automation boundaries and stilde by maximizing net output
    % start maximization at the point where automation first hits:
    startauto = [rno.sgrid(indmin) rno.x(indmin) rno.x(indmin)];
    % net output maximization:
    %[eq,~,ef] = fmincon(@(eq) -output(sbar,eq,p),startauto,...
    %    [0 1 -1],0,[],[],[0 0 0],[1 1 1],[],options);
    % including the constraint that output must be strictly positive, to
    % exclude local minimum at zero (may be necessary for some parameter
    % constellations, uncomment if needed):
    [eq,~,ef] = fmincon(@(eq) -output(sbar,eq,p),startauto,...
     [0 1 -1],0,[],[],[sbar 0 0],[1 1 1],@(eq) outputcon(sbar,eq,p),options);
    
    % compute labor assignment and wages
    rlab = laborassignment(sbar,eq(1),eq(2),eq(3),p);
    
    % collect results
    r = rlab;
    r.stilde = eq(1);
    r.xlow = eq(2);
    r.xup = eq(3);
    r.qzero = qzero;
    r.exitflag = (ef==1 | ef==2);
    r.exitflag1 = "interior automation";
    
elseif (compadsbar == 1) && (compadmostskilled == 1) %interior or low skill
    
    % compute assignment, wages, output assuming low-skill automation
    xup = automationup(sbar,sbar,0,p);
    rlab = laborassignment(sbar,sbar,0,xup,p);
    
    % check whether skill sbar is more expensive than capital in task 0
    if rlab.w(1)*exp(logcapprod(0,p))/exp(loglaborprod(sbar,0,p)) >= 1 % low-skill automation
        
        r = rlab;
        r.stilde = sbar;
        r.xlow = 0;
        r.xup = xup;
        r.qzero = qzero;
        r.exitflag = 1;
        r.exitflag1 = "low-skill automation";
        
    else % interior automation
        
        % compute automation boundaries and stilde by maximizing net output
        % start maximization at the point where automation first hits:
        startauto = [rno.sgrid(indmin) rno.x(indmin) rno.x(indmin)];
        % net output maximization:
        %[eq,~,ef] = fmincon(@(eq) -output(sbar,eq,p),startauto,...
        %    [0 1 -1],0,[],[],[0 0 0],[1 1 1],[],options);
        % including the constraint that output must be strictly positive, to
        % exclude local minimum at zero (may be necessary for some parameter
        % constellations, uncomment if needed):
        [eq,~,ef] = fmincon(@(eq) -output(sbar,eq,p),startauto,...
         [0 1 -1],0,[],[],[sbar 0 0],[1 1 1],@(eq) outputcon(sbar,eq,p),options);
    
        % compute labor assignment and wages
        rlab = laborassignment(sbar,eq(1),eq(2),eq(3),p);
        
        % collect results
        r = rlab;
        r.stilde = eq(1);
        r.xlow = eq(2);
        r.xup = eq(3);
        r.qzero = qzero;
        r.exitflag = (ef==1|ef==2);
        r.exitflag1 = "interior automation";
        
    end
    
elseif (compadsbar == 0) && (compadmostskilled == 0) %interior or high skill
    
    % compute assignment, wages, output assuming high-skill automation
    xlow = automationlow(sbar,1,1,p);
    rlab = laborassignment(sbar,1,xlow,1,p);
    
    % check whether skill 1 is more expensive than capital in task 1
    if rlab.w(end)*exp(logcapprod(1,p))/exp(loglaborprod(1,1,p)) >= 1 % high-skill automation
        
        r = rlab;
        r.stilde = 1;
        r.xlow = xlow;
        r.xup = 1;
        r.qzero = qzero;
        r.exitflag = 1;
        r.exitflag1 = "high-skill automation";
        
    else % interior automation
        
        % compute automation boundaries and stilde by maximizing net output
        % start maximization at the point where automation first hits:
        startauto = [rno.sgrid(indmin) rno.x(indmin) rno.x(indmin)];
        % net output maximization:
        %[eq,~,ef] = fmincon(@(eq) -output(sbar,eq,p),startauto,...
        %    [0 1 -1],0,[],[],[0 0 0],[1 1 1],[],options);
        % including the constraint that output must be strictly positive, to
        % exclude local minimum at zero (may be necessary for some parameter
        % constellations, uncomment if needed):
        [eq,~,ef] = fmincon(@(eq) -output(sbar,eq,p),startauto,...
         [0 1 -1],0,[],[],[sbar 0 0],[1 1 1],@(eq) outputcon(sbar,eq,p),options);
        
        % compute labor assignment and wages
        rlab = laborassignment(sbar,eq(1),eq(2),eq(3),p);
        
        % collect results
        r = rlab;
        r.stilde = eq(1);
        r.xlow = eq(2);
        r.xup = eq(3);
        r.qzero = qzero;
        r.exitflag = (ef==1|ef==2);
        r.exitflag1 = "interior automation";
        
    end
    
else %incompatible with comparative advantage across workers
    
    r.exitflag = 0;
    r.exitflag1 = "Comparative advantage assumption across workers is violated";
    return
    
end

% provide additional information about possible automation patterns
if compadleastskilled==0 && compadmostskilled==1
    r.exitflag2 = "w/o minimum wage: automation is always interior";
elseif compadleastskilled==1 && compadmostskilled ==1
    r.exitflag2 = "w/o minimum wage: automation transitions from interior to low skill";
elseif compadleastskilled==0 && compadmostskilled==0
    r.exitflag2 = "w/o minimum wage: automation transitions from interior to high skill";
end

end